clear; clc; close all;

%% 固定参数
a0 = -2;
a1 = 0.83;
a2 = 0.14;
b0 = 1.3;
b1 = 0.1;
b2 = 1;
d1 = 0;
d2 = 0;

% y0, z0 固定
y_fixed = 0.1;
z_fixed = 0.1;

x0_min = -60;  
x0_max = 60; 
num_x0 = 1201;   % 例如间隔为 1
x0_values = linspace(x0_min, x0_max, num_x0);

%% ========== 迭代设置 ==========
N_total = 1000;    % 总迭代步数
N_trans = 200;    % 舍弃瞬态步数
sample_interval = 10;  % 分叉图采样间隔

%% ========== 存储结果的数组 ==========
LE_result = zeros(num_x0, 3);  % 每个 x0 下的 (LE1, LE2, LE3)
x0_bif = [];  % 分叉图横坐标
x_bif  = [];  % 分叉图纵坐标 (记录 x_n)

%% ========== 扫描 x0 ==========
for i = 1:num_x0
    x_init = x0_values(i);
    
    % 初始状态
    x = x_init;
    y = y_fixed;
    z = z_fixed;
    
    % Lyapunov 指数算法初始
    Q = eye(3);
    sum_log = zeros(1,3);
    
    % 迭代计算
    for n = 1:N_total
        % 一步离散映射
        x_new = x + y;
        y_new = sin(a0 * y^2 * sin(a1 + a2*x)) + b2*sin(z) + d1;
        z_new = b0*y + b1*z + d2;
        
        % 计算当前状态下 Jacobian
        %  J = [1, 1, 0;
        %       M1, M2, b2*cos(z);
        %       0,  b0, b1]
        theta = a1 + a2*x;
        sin_theta = sin(theta);
        cos_theta = cos(theta);
        arg = a0 * y^2 * sin_theta;
        cos_arg = cos(arg);
        
        M1 = a0 * a2 * y^2 * cos_theta * cos_arg;
        M2 = 2 * a0 * y * cos_arg * sin_theta;
        
        J = [1, 1, 0;
             M1, M2, b2*cos(z);
             0,  b0, b1];
        
        % QR 正交化
        Q = J * Q;
        [Q, R] = qr(Q);
        % 累计 log|Rii|
        if n > N_trans
            diagR = diag(R);
            sum_log = sum_log + log(abs(diagR))';
        end
        
        % 状态更新
        x = x_new;
        y = y_new;
        z = z_new;
        
        % 分叉图数据记录：舍弃瞬态后每 sample_interval 记录一次
        if (n > N_trans) && (mod(n, sample_interval)==0)
            x_bif(end+1) = x;       % 记录当前 x
            x0_bif(end+1) = x_init; % 横坐标是这个 x0
        end
    end
    
    % 计算 LE
    num_effective = N_total - N_trans;
    LE = sum_log / num_effective;
    LE = sort(LE, 'descend');  % LE1 >= LE2 >= LE3
    LE_result(i,:) = LE;
    
    % 显示进度
    fprintf('Progress: %.1f%% \n', (i/num_x0)*100);
end

%% ========== 绘图 ==========
figure;

% 上子图：Lyapunov 指数 vs x0
subplot(2,1,1);
hold on; grid on;
plot(x0_values, LE_result(:,1), 'r-', 'LineWidth', 1, 'DisplayName','LE_1');
plot(x0_values, LE_result(:,2), 'b-', 'LineWidth', 1, 'DisplayName','LE_2');
plot(x0_values, LE_result(:,3), 'g-', 'LineWidth', 1, 'DisplayName','LE_3');
plot(x0_values, zeros(size(x0_values)), 'k--', 'LineWidth', 1, 'DisplayName','0');
xlabel('x_0','FontSize',12);
ylabel('LE','FontSize',12);
title('x_0 李指数图','FontSize',14);
legend('Location','best');

% 下子图：分叉图 ( x_n vs. x_0 )
subplot(2,1,2);
plot(x0_bif, x_bif, 'r.', 'MarkerSize',1);
grid on;
xlabel('x_0','FontSize',12);
ylabel('x_n','FontSize',12);
title('x_n - x_0 分岔图','FontSize',14);
